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SQIF  Array  Dynamics 


1.  Abstract 

We  develop  an  analytic  perturbation  solution  for  single-cell  SQUIDs  and  for  multi-cell  rows 

with  random  cell  sizes  (constrained  by  having  a  mirror  symmetry  about  the  center  junction).  The 
response  of  a  single  SQUID  to  an  external  RF  signal  is  determined  in  terms  of  correlation 
functions  between  the  phase  of  the  inline  effective  junction  and  the  circulating  current.  This 
yields  an  RF  phase  shift  that  can  be  calculated  as  a  function  of  the  circuit  parameters  and  the 
signal  frequency.  Summing  the  phase-shifted  RF  voltages  of  many  random-cell-size  SQUIDs  in 
series  yields  the  SQIF  RF  response. 

2.  Characterization  of  SQIF’s  as  RF  Detectors 

SQIF  arrays  were  originally  developed  to  detect  extremely  small  static  and  quasistatic 
magnetic  fields.  The  paper  by  Schopohl  et.  al.1  showed  how  it  was  possible  to  approximate  the 
RF  response  of  a  SQIF  array  to  a  low  frequency  RF  signal  by  using  the  transfer  function  (DC 
voltage  vs  external  magnetic  field)  in  the  narrow  region  near  zero  magnetic  field,  where  the 
transfer  function  is  essentially  linear.  In  a  SQIF  array  the  interference  of  the  voltages  of  the 
random-area-SQUID  loops  in  the  array  leads  to  the  steep  slope  of  the  transfer  function  near  the 
origin,  which  should  be  the  only  place  where  the  voltages  interfere  constructively. 

To  characterize  the  merit  of  a  SQIF  array  as  an  RF  detector,  it  is  necessary  to  take  into 
account  both  the  Johnson-Nyquist  noise  of  the  array  junctions  and  the  appearance  of  intrinsic 
nonlinear  effects  which  determine  the  range  of  signal  amplitudes  in  which  the  transfer  function 
can  be  considered  linear.  The  widely  accepted  figure  of  merit  for  RF  detectors  is  the  ‘spur  free 
dynamic  range’  (SFDR),  which  is  a  measure  of  the  range  of  detector  linearity,  and  is  given  by  the 
range  (in  dB)  where  the  detector  response  (v2)  for  two  signals  that  are  close  in  frequency  is  above 
the  Johnson-Nyquist  noise  threshold,  while  the  slightly  frequency-shifted  response  due  to  cubic 
mixing  (at  a)=20)i-a)2  or  co=2co2-tOi  )  remains  below  the  noise  threshold. 

The  approach  that  we  have  taken  is  to  fully  determine  the  RF  response  of  each  row  of  a  SQIF 
array  to  the  magnetic  component  of  an  RF  field,  and  sum  the  voltages  to  obtain  the  total  RF 
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response.  This  means  that  both  the  amplitude  and  RF  phase  of  the  response  of  each  row  at  the 
signal  frequency  are  determined,  in  contrast  to  the  approach  in  [1]  where  the  phase  shifts  of  each 
row  are  taken  to  be  the  phase  shifts  at  DC  (i.e.  those  due  to  the  static  external  flux)  which  are 
then  assumed  to  be  applicable  to  low  frequency  signals.  The  approach  in  [1]  essentially  treats  the 
RF  signal  as  an  adiabatically  changing  flux  bias. 

To  determine  the  SFDR  of  an  array  we  start  with  the  response  of  a  single  row,  and  with  the 
simplest  case  of  a  single  SQUID  or  a  single  row  SQIF  with  zero  inductance  that  can  be  treated  as 
an  equivalent  effective  single  SQUID1.  The  linear  response  of  the  effective  SQUID  when  the 
input  is  a  pair  of  small  amplitude  signals  (amplitudes  {arfl,aif2}  at  frequencies  {(0Ilv(0lf2 }  )  can 

be  solved,  with  solutions  in  the  form  of  a  sum  of  sinusoids  at  the  signal  frequencies  o)rf,i  with 
amplitudes  A;  and  phase  shifts  HQ ,  .  For  high  Tc  junctions  /3r  ~  0  ,  and  this  allows  good 

approximate  solutions  to  be  obtained  for  the  unperturbed  SQUID  voltage  £  and  induced  flux 
rjDC  =PJdrc  >  where  the  DC  subscript  refers  to  the  response  to  <f)DC ,  the  DC  component  of  the 
magnetic  flux  threading  the  SQUID  loop  (the  /3C  =0  solution  will  be  discussed  further  below).  If 
the  linear  response  (S^RF)  is  assumed  to  be  a  valid  approximation  for  small  amplitude  RF  signals, 
then  it  can  be  determined  in  terms  of  the  correlation  functions 


sin(Qc)sin( 


;„/  *Pdc  Vdc 


)j  and  (cos^dc)cos^dc+2,1dc^ 


as 


%KF  =  X  A  sin  (<*V,T  -  ’*'*.,)=  X 


2  fa  ^ 

Urf,i 


V  2  J 


sin(CDC)sin(^A^) 


®Q+  (cos(CDC)cos(^U^^) 


,sin(ffl7/iT  —  HQ,.), 


where  the  RF  phase  shifts  are  given  by 


HQ,  ~  sin  1 


qv,,-  /  Jaiff(c°s(CDc)cos(^DC  t 71 DC ) 


(l) 


(2) 


(3) 


It  can  be  seen  that  as  long  as  the  (cos cos)  correlation  function  is  non-zero,  the  RF  phase  shift 
HU  will  approach  zero  as  the  signal  frequency  corj  I  goes  to  zero.  A  useful  corollary  to  this 
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behavior  is  that  as  the  (cos cos)  correlation  function  becomes  small  and  of  the  same  order  as 
coif ,  the  RF  phase  shift  will  depend  sensitively  on  [3,  or  <f)DC  for  fixed  (Orf , ,  or  conversely  it 
will  be  a  steep  function  of  cor/  i  when  f3L  and  0DC  are  fixed.  It  will  be  shown  below  that  the 
correlation  function  (cos  cos)  becomes  small  when  either  /3L  or  (j)DC  is  small. 

3.  Derivation  of  the  Hamiltonians  for  Mirror  Symmetric  Row  Arrays 

The  row  arrays  that  we  shall  treat  here  can  be  considered  as  a  subset  of  the  completely  random 
SQIF  array  row  described  by  Oppenlander  et  al.1  in  that  they  are  random  but  have  mirror 
symmetry  about  a  central  junction.  The  Hamiltonians  and  equations  of  motion  for  mirror 
symmetric  row  arrays  are  described  by  Ncell  + 1  junction  phases,  along  with  Ncell  constraint 

equations,  and  Ncell  circulating  currents  (or  their  equivalent  fluxes),  for  a  total  of  NceU  + 1 

independent  coordinates.  The  Ncell  constraints  may  be  included  in  the  Lagrangian  along  with 

N cell  Lagrange  multipliers  before  forming  the  Hamiltonian.  However,  a  simpler  method  is  to 

transform  to  a  set  of  NceU  + 1  coordinates  that  include  the  constraints,  which  then  yield  one 

independent  junction  phase  and  Ncell  independent  circulating  currents.  We  illustrate  this  with  an 

NCeu  =4  array  shown  below,  where  the  middle  junction  is  allowed  to  have  a  critical  current  that 

differs  from  the  remaining  junctions  by  a  factor  of  p. 


Figure  1.  Four-Cell  Array 
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We  first  transform  the  external  fluxes  (0i>  02  >  03’  04 }  to  the  set  {0fl,  0/>-  0C>  0  s !  ,  where 
0S  =  (0,  +  02  +  0:,  +  04 )  /  4  would  be  symmetric  in  all  four  fluxes  (if  the  fluxes  were  equal),  while 
the  remaining  iVceH  - 1  coordinates  are  the  partly  antisymmetric  combinations 

t  =  (01  +  02  -  -  04  )  /  4 

=  (01  -  02  -  03  +  04  )  /  4 

(4) 

<Pc  =  (01  —  02  +  03  -  04)  /  4- 

The  iVceH  induced  fluxes  fl/(  =PLMicirCrk  are  transformed  in  an  identical  way.  By  writing 
=  £  +  a;.(0a  +ria)  +  bj((j)b  +Tjb)  +  Cj((l)h  +  ric)  +  Sj((f)s  +77J,  the  constraint  equations  for  the  four 
cells,  along  with  the  demand  that  the  coefficients  of  the  four  cross  terms  C0,  in  the  kinetic  energy 
be  zero,  can  be  solved  to  yield  a  solution  for  the  coefficients  [  a; ,  /?, ,  c(. ,  y  (  .  The  transformation 
of  the  original  junction  phases  are  then  given  by 


4  =  C+2 

4  =  C+ 
4  =  C  - 
&=C  + 


4  =  C+2 


(P+i) 
(P  +  4) 
(P-2) 
(P  +  4) 
6 

(P  +  4) 
(P-2) 
(P  +  4) 
(P  +  1) 
(P  +  4) 


2 

(<t>a  +PJ  +  2(PS  +  T  )  +  (0C  +  T),  ) 

(0.  +  na ) + (0,  +  ps )  -  (0ft  +  pj- 


(0u+pj  + 

(0«  +  Va  )  -  (0S  +  P  y  )  +  (0/,  +  flfc  )  “ 


2 

(P  +  4) 
(P  +  2) 
(P  +  4) 


(0fl+Pj-2(0,+Ps)  + 


2 

(P  +  4) 


(0f +  *U 


(0,  +  P ) 

(0<  +  P) 

(0c  +  Pc) 


(5) 


It  can  be  seen  that  the  phase  of  the  middle  junction  will  have  contributions  from  the  fluxes 

that  are  antisymmetric  (if  all  cell  fluxes  are  equal)  about  the  center  conductor  in  addition  to  that 
from  the  effective  junction  phase  C  • 

When  p  =1  the  effective  phase  C  may  be  obtained  by  summing  all  of  the  junction  phases  and 
dividing  by  five.  When  p  ^  1 ,  the  expression  for  C  will  also  include  flux  terms  multiplied  by  a 
factor  proportional  to  (p  —  1) .  The  expression  for  the  junction  potential  energy  may  be  written 
as 
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VJJ  /  EJ  =  -%  cos(4 )  -  (p  - 1)  cos(4 ) 

i= 1 

=  -2  cos[C  +  2^(0a  +  pa)  +  +  nc )]  cos[2(0s  +  ns )] 

-  2  cos[£  +  (0fl  +  rja  )-££(&  +  Uc )]  cos[0,  +  r\s )  -  (0/,  +  Vb )] 

-  p  cos[C  -  ^4  (0fl  +  p  J  +  ^4  (0e  +  pc )]. 


For  the  case  of  a  two-cell  SQUID,  the  transformed  coordinates  and  junction  potential  function 
are  given  by 


£  =  C  +  (0,  + )  +  — 7T7  (0a  +  Va ) 
(P  +  2) 


&=c- 


(P  +  2) 


(<Pa+rla) 


$3  =  C-  (0.v  +  )  +  —777  (0fl  +  ), 

(P  +  2) 


(7) 


and 


Vjj  / Ej  =  -^cos(^)-(p-\)cos^ ;2) 

i=l 

=  -2  cos[f  +  (0a  +  Tla  )]  COS[(0,  +  ps )] 

-  pCOSiC- ^2  (0a+ ^il¬ 


ls) 


For  the  case  of  a  one-cell  SQUID,  the  transformed  coordinates  and  junction  potential  energy 
are  given  by 

4  =  C+(0+p)/2 

4=C-(0+77)/2,  (9) 


and 


Vjj/Ej=-%c°s(Zi  )  =  -2  cos(£)cos[(0  +  p)  /  2], 

i=l 


(to) 
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4.  fic  =  0  Solution  for  a  Single  SQUID 


Using  the  coordinates  given  above  in  Eq.  9,  the  kinetic  energy  term  in  the  Hamiltonian  may  be 
written  as 


T 


C 

2 


fo  V 

_ 0_ 

v  2  K  ; 


{2C2+t);/2}, 


(11) 


and  in  dimensionless  units  as 
T  M(  M 

—  =  ^fC  +  ^r7lr 

E,  2  2 


with 


M(=2fic,  Mn=ficl  2 

d  Pjcirc  ’ 

Ej  =Ic®q/2k, 


where 


I3C=C 


(  2n^ 

V°o  j 


2  /-»  2/r 

IcRn>  Pl  —  L(2kIc /O0),  t  —coj,  C0c  —  ^  IcRn . 


The  inductive  potential  energy  may  be  written  as 


K  =L/1  /  2  =Ey/3Lil  /  2  =£,  T]2 


with 


CO, 


PcPl 


’  Eire  l circ  ^  ^ c  * 


The  total  Hamiltonian,  including  the  forcing  current  potential,  is  then  given  by 


H=T  +  VJJ  +  Vll-EJib£, 

which  yields  the  equations  of  motion  in  terms  of  the  transformed  coordinates  {£,  t?}  as 


(12) 


(13) 


(14) 

(15) 

(16) 

(17) 


Pc  C’TT  +C>T  +sin(C)cos[(0  +  77)  /  2]  =  ib 

Pc  +t?,T  +2 cos(C)sin[(0  +  77)  /  2]  +  =  0  .  (is) 
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For  over-damped  junctions  we  set  fi(  —0  and  rewrite  C,  as  C  =  0  +  K  /  2 .  The  equations  of 
motion  for  0  can  then  be  written  in  the  form  used  by  Swift  et  al.1  as 


=  1  -  cos(0)cos[(0  +  77)  /  2]  /  4  =  1  +  a  cos (0) 

77, T  +(2  /  [3L  )Tj  =  2  sin(0)sin[(0  +  77)  /  2],  (19) 

where 

fl  =  -cos[(0  +  77)  /  2]  /  4  and  z  =zib.  (20) 

Swift,  Strogatz,  and  Wiesenfeld2  have  shown  that  for  a  system  of  single  junctions  connected  in 
series,  a  perturbation  theory,  based  on  the  exact  solution  for  a  single  junction,  can  be  developed 
to  determine  the  DC  voltage  for  the  coupled  junction  system.  The  set  of  relative  phases  for  the 
junctions  is  not  determined  by  the  theory,  but  there  are  optimal  sets  where  the  stability  of  the 
array  is  maximized.  For  the  SQUID  problem,  we  have  extended  this  theory  to  determine,  in 
addition  to  the  perturbed  frequency  (yDC[4,0DC,/3L],  the  time  dependence  of  the  effective 

junction  phase  C,DC  =^0(coDCt),  and  the  induced  flux  riDC  =P,icirc/)C  (or  equivalently  a  description 

in  terms  of  the  two  junction  phases  and  q2 )  as  functions  of  ih,  0  and  jiL . 

To  develop  a  perturbation  expansion  for  0  and  V ,  we  rewrite  the  0  equation  of  motion  as 


0,Y  =  l +  a()  cos(0)  +  (a-  a0 ) cos(0) ,  with  a0  =  -cos[(0  +  (77)) / 2] / 4,  (21) 

where 

1  r 

(77)  =  — [  ri(coz)dT,  T  =  2k/co,  77  =  (77)  +  <577.  (22) 

J  JO 

Treating  the  last  term  as  a  perturbation,  the  unperturbed  equation  of  motion  for  0  may  be 
solved  exactly  yielding 

0o(T)=2tan-‘(Vtan((7)o/2)),  (23^ 


where 


<Pr 


\j  1  £7()  z  -\j  1  ct q  /.,  T  (t)i  j  T  and  ^0 


1  +  an 


(24) 
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To  obtain  the  first-order  correction  to  the  frequency,  we  use  the  0  equation  of  motion  to 
write 


co  _d(p_  dcp  cW  ^/T-i 


4  dz  clO  dz  1  +  a0  cos(6) 


[l  +  a0  cos (6)  +  (a-  ao)cos(0)] 


:  ~  Gn  +  - 


1  +  a0  cos (6) 
Using  the  trigonometric  identity1 


(a  -  a0 )  cos(d)  =  ( co0  +  8co)/ib. 


(25) 


2 _ =  1  ~  a” c°s(y»)  and  cos(8.)  =  cos(«2-fl' 


1  +  a0  cos (64 ) 


1 


1  -  a0  cos(<p0 ) 


(26) 


the  expression  for  the  first-order  frequency  correction  may  then  be  written  as 

1  1  f 


8co  =4 


27r 


J0  (cos ((p)-a0)(a((p)-a0)dcp. 


(27) 


To  proceed  further  it  is  necessary  to  solve  the  equation  of  motion  for  q  to  lowest  order,  which  is 
obtained  by  expanding  the  sin[(0  +  77)  /  2]  term  in  the  equation  of  motion  to  first  order  in  77 ,  with 
solutions  given  the  label  p0 

77,  r  +(2  /  >8^  )?7  =  2  sin(0o )  sin[T  (^>  +  rj)]  =  2sin(0o){sin({f)(l-{772  +  ...)  +  cosf^Kfp-...} 
rio’r+(2/ «2sin(0o){sin(^)  +  cos({0H?7o}.  (28) 

Using  the  exact  zero-order  solution  for  6>„,  sin(0o)  and  cos(0o)  may  be  evaluated  as 


l-«o  sin(ft»0T) 


sin(e,)  = - ,  cos(e.)  -  cos(ffl» T)‘a» 


1  -  a0  cos(ft)0  t)  u'  1  -  a0  cos(ft)0  t) 

This  allows  the  re-expression  of  the  lowest  order  77,,  equation  of  motion  as 


(29) 


770  T  +  /(f)[sin(0  /  2)  +  770  cos(0  /  2)]  +  2rj0  /  pL  =  0, 


where 


/(T) 


.  sin(W0T)/s,n 

1  -  COS(ft>0T) 


,  /sin  =-2^. 


(30) 


(31) 


As  a  first-order  ordinary  differential  equation,  the  77,,  equation  of  motion  may  be  solved  exactly 
as 
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r£(j)  =  -e 


-Jo'(2/ft+cos(0/2)/(T'))</T'  fl 


Jorsin(0/2)/(T')e 


f  (2/i3i+cos(0/2)/(T',))rfT' 
'Wj  Jo 


dr'. 


The  integral  in  the  two  exponential  factors  may  be  carried  out  to  yield 


f  f(z ')clT '  - a°  cos(®oT)] - / 2 } 


a0®o 


allowing  the  solution  to  be  written  as 


__  -2t/j8l  (V*-“o  cos(0/2)/aowo)ln[l-aocos(®oT)]  . 

1 10  e  z 


g2z'l[lL ^(-yi-ao  cos(^/2)/a0co0)ln|l-«0cos((00r')] 


/sin  sin(0  /  2)sin(co0Tr)  ,  .  , 


(l-a0  cos((y0r')) 


r/r'. 


Due  to  the  fact  that 


(32) 


(33) 


(34) 


■\jl-al  cos (0 12)1  a0(O0 )  =-l, 


(35) 


the  denominator  in  the  integrand  is  canceled  out  and  the  integration  can  be  carried  out  exactly, 
yielding 


ty,/) = ' 


—  sin(0/2) 


where 


c  =  2/(co0Pl),  k  =  - 


cos(«0t)  -  c  sin((n0T) 
(1  +  c2  )(1  -  aQ  cos(ft)0T)) 


2  sin(0/2) 

4  1  +  c2 


=  ?c 


cos(«0t)  -  c  sin(tu0T) 
l-a0  cos(£00t) 


The  first-order  correction  to  the  rotation  frequency 


(36) 


(37) 


8(0  =i. 


1  co0  r2rc 


r  in 

J o  (cos (ffl0T) -  a0)[a((p0 ) -  aQ (q o0 )] d z, 


(38) 


may  now  be  carried  out  by  expanding  a((p0)  in  terms  of  powers  of  /?,, .  The  integration  required 
after  expanding 

a(<Po)  =  -  cos  [j(0  + 1?0/))]  /  4  (39) 

to  any  desired  order  in  q0 ,  may  be  carried  out  in  closed  form.  To  second  order  in  Vo  this  yields 
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Making  the  replacement 

co0  — >  coDC  =  co0  +  8(0  (43) 

in  the  expressions  for  1} 0  and  d0  then  yields  the  first  order  unperturbed  solutions  (jDC  and  dDC . 
The  unperturbed  flux  r\DC  can  also  be  expressed  in  terms  of  Co  as 

Inc  =  -  ^  +  c\l-a0) sin(Co  -(n/2  +  yf)) 

(44) 

=  p  +  q  sin(£0  -  (7T  /  2  +  yr)), 


where  the  phase  shift  V  is  given  by 


Vr  =  sin'  (i'J 3L/2)/ 


(45) 


In  this  form,  Vdc  has  both  a  constant  and  an  oscillating  component,  and  the  result  can  be  seen  to 
make  sense  physically  by  considering  the  supercurrent  Ic  sin(C)  going  through  the  effective  C 
junction.  This  current  amplitude  will  be  a  maximum  when  £  =  ±^r  /  2,  at  which  point  r lDc  will 
be  at  a  minimum  (for  y/  =0 )  and  vice  versa.  As  the  SQUID  inductance  increases  from  zero,  the 
oscillating  part  of  Vdc  will  be  shifted  in  phase  by  an  amount  that  is  a  function  of  the  flux  ih[3, . 


Below  is  a  plot  showing  TjDC  (black)  and  rjDC  obtained  by  a  numerical  solution  (green)  of  the 
coupled  ODEs  for  ib  =1.05,  (pDC  =0.2/r  and  fiL  =1.0,  where  the  motion  is  markedly  non- 
sinusoidal.  It  can  be  seen  that  the  value  obtained  for  the  frequency  is  fairly  accurate,  while  the 
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analytic  expression  for  the  amplitude  of  Vdc  differs  from  the  amplitude  of  the  numerical 
solution  by  about  one  to  two  percent. 


n 


Figure  2.  SQUID-Induced  Flux  /]( t)  vs  t 


5.  Linear  Response  of  a  SQUID  to  an  RF  Signal  for  the  (5C  =0  Case 

To  determine  the  linear  response  of  a  SQUID  to  an  rf  signal,  we  write  the  equations  of 
motion  for  the  effective  phase  C 

A  C’rr  +sin[C]cos[(0  +  77)  /  2]  =  4  (46) 

in  terms  of  the  unperturbed  solution  { £DC ,  /)/)C}  and  denote  the  linear  response  to  the  RF  signal 
as  <54 ,  where  the  RF  signal  is  given  by  8(j)  =  A  sm(corfz).  Rewriting  the  equations  of  motion  with 

C  =  Cdc  +  <5C f  and  0  =  0 DC  +  A  sin(m^T)  (47) 

then  yields 

0C,/,r  +  cos[4DC  ]  cos  [(<pDC  +  r\DC )  /  2]^  -  sin[4DC  ]  sin[(0oc  +  r\DC )  /  2]  f  sin(mr/r)  =  0,  (48) 

where  we  have  subtracted  out  the  DC  solution  and  kept  only  linear  terms  in  8C,  and  A.  Assuming 
that  the  frequency  of  the  RF  signal  is  much  less  than  (0DC  ,  it  is  reasonable  to  average  over  the 
rapid  oscillations  of  and  t]dc  which  yields 

0C,/>r  +  (cos[CDC]cos[(0DC  +r\DC)l  2])8Crf  -  (sin[£DC]sin[(0DC  +  rjDC)  /  2])f  sin^r)  =  0.  (49) 
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The  solution  for  the  linear  response  SCIU  can  be  written  in  terms  of  the  two  correlation  functions 
as 


=  f (sin^oc)sin( 


hoc 


s(Cdc  )  c°s(  ^DC  +2  7100  )j  sin(ft>^T)  -  cos(ft)^T) 

cos(CDC)cos  (^77^)V  +  < 


(50) 


This  may  be  re-expressed  in  the  simpler  form  of  a  single  sinusoid  as 


<5£/  =  7^in(£DC)sm(^DC*??DC 


sin(<o „) 


ft)‘+(cos(fDC)c°s( 


fipc  hoc 


(51) 


where  the  RF  phase  shifts  is  given  by 


Tr/  =  sin1 


7  J®^+  lcos(^DC)cos(^DC  j '1dc 


(52) 


V  '7 

It  can  be  easily  seen  that  'F,/  will  approach  zero  as  the  signal  frequency  co,f  goes  to  zero.  A 
useful  corollary  to  this  behavior  is  that  as  the  (cos  cos)  correlation  function  becomes  small  and 
of  the  same  order  as  (Olf ,  the  RF  phase  shift  will  depend  sensitively  on  [5L  or  (pDC  for  fixed  corf , 
or  conversely  it  will  be  a  steep  function  of  co^  when  [3L  and  <f)DC  are  fixed. 


By  using  the  expressions 


cos(Cdc)=-A  a;x;Sin(0J,x'T),  sin(Cgc)  =  COS(ft)gcT)  a°c  , 

l-aDCcos(eoDCz)  1  -  aDC  cos(coDC  z)  (53) 

the  correlation  functions  can  expanded  to  any  order  in  (t?DC)"  and  integrated  term  by  term.  Thus 
to  second  order  in  TjDC  we  have 


•  y  CO  nr  C^k^(0dc  ^ 

(c°s[Cdc]cos[(0dc  +  t]Dc)  /  2]y  —  J0  cos[£dc(t)] cos [{(/)dc  +  ^?£>c(^))  !  2]^T 


(54) 


CO 


DC 


sin(0s.  /  2) _ (1  -I-  aDC 

- IV  C 


& nr (1  anc^ 


v)  ,  cos (<ps/2)  ^_ 

“I  IV  c 


{1  bDC+aDC\  1- 


■  aDC{ 1  +  a 


DC 


1b~l 

2  UDC 


)]} 


(1  cinc)ciDr(\-\- 

anc ) 


and 


(55) 
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(sin[Cz5c]sin[(0flc  +  Vdc)  1 2l)  =  f  *  sin[^oc(T)]sin[((/)oc  +  rjDC(z))  /  2]dt 

in  J0 

-  cos (</>s  /  2)sin(0s  /  2)/,  +  (cos (0s  /  2)2  +  cos(05  /  2)  - l)^/2  -  aflC[cos(0s  /  2)  -  }sin (0s  /  2)2]^/3 

2  2  2 

-sin(0s  /2)cos(f  /2)^-(l-c2)/4  +aDC sin(0s  /2)cos(0s  /2)^-/s  -sin(0s  /2)cos(0s  /2)^-c2I6. 

o  o  o 

The  expressions  for  the  integrals  7,  through  /6  are  given  in  appendix  A. 

Plots  of  the  two  correlation  functions  ^sin(£oc)sin(^oc  +  ^DC  )j  and  (cas(CDC)cos(^'>c  * 11,x  )^j, 

expanded  to  second  order  in  /7dc  ,  are  plotted  below  for  ib  =1.05  and  4  =4  (n.b.  vertical  scales 
differ). 


(cos(£Dc)cos((0Dc+/7dc)/2))  (4=1.05)  <cos(£nr)cos((<Ar>r+f7nr)/2)>  (4=4.0) 


vs.  0DCand  f$L  with  ib  —1.05  vs.  0/x-  and  with  ib  —4.0 
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It  can  be  seen  that  for  small  values  of  or  PL  for  any  value  of  4 ,  or  for  large  PL  and  large  4 , 
that  (cos cos)  will  be  small,  allowing  the  RF  phase  shift  to  be  sensitive  to  small  changes  in 
frequency.  We  have  plotted  the  RF  phase  shift  below  for  =  5.0,  for  signal  frequencies  from 
(0rf  =0  to  0Jrf  =0.01  ft)  r .  It  can  be  seen  that  vFr/  is  a  very  steep  function  of  corf  for  small  values 

of  (f)DC  for  the  chosen  SQUID  parameters  ( ib  =4  for  both  plots). 


Figure  7. 

Two  Cell  SQUID  Phase  vs  (pl)C  and  (0^ 


Figure  8. 

Two  Cell  SQUID  Phase  vs  (pl)C  and  C0rf 


14 


6.  /3C=0  Solution  for  a  Two-Cell  SQUID 


For  external  magnetic  fluxes  near  (f)a  ~  0,  (j)s  ~  0 ,  the  equations  of  motion  for  the  two-cell 
SQUID  can  be  solved  using  the  perturbation  method  developed  above  for  the  single  cell  SQUID. 
From  the  numerical  solution  for  the  two-cell  case,  it  can  be  seen,  however,  that  using  the 
approximate  solution  of  Eq.  23  for  the  effective  £  junction  (which  consists  of  a  string  of  single 

pulses),  will  not  be  valid  when  the  t]a  and  r\s  amplitudes  become  large  enough  to  markedly  affect 
the  C  junction,  as  shown  in  the  example  below. 


^DC  Voltage  (4=1.05,  <pa=3.U,  tps=  1.04)  r/ax  (4=1.05,  0«=3.14,  0S=1.O4) 

bDC,T  f}aDC 


Two-Cell  SQUID  Voltage  d£DC(z)  /  dz  vs.  T  Two-Cell  SQUID  Induced  Flux  T]^  (z)  vs.  Z 

r]sX  (4=1-05,  0a=3.14,  <A,=  1 .04( 


TJsdc 


Two-Cell  SQUID  Induced  Flux  (Z)  vs.  Z 
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To  cover  the  entire  range  for  0„  and  <ps  it  will  be  necessary  to  include  the  time-dependent 
tj  and  rjs  solutions  in  determining  C(T)  and  perhaps  an  iteration  of  the  coupled  solutions.  In 
what  follows,  we  develop  the  solutions  for  the  case  where  0n  and  0t  are  near  zero. 

Using  the  coordinates  given  above  in  Eq.  7,  the  kinetic  energy  term  in  the  Hamiltonian  may 
be  written  as 


2 


f  <j>  ; 

2( 

0 

k2k  y 

1 

2  ,  ,  2p  2 

la 


(2  +  p)C  +2t)s  + 


2  +  P  J  , 


and  in  dimensionless  units  as 

T  M(  Ms  2  M 
-  =  -r^  +  -fn-,T  +  -rri[ 


la,  T 


with 


Mt.  =  (2  +  p)j3c,  M  =  2ft.,  M»=^-Pc 
V,=PJ„  Va  =  PJa>  P 

is  =  (it+i2)/2,  4  =  0;- i2)/ 2, 

Pls=Pl+Pm >  pLa=PL-pM- 


(56) 


(57) 


(58) 


The  inductive  potential  energy  may  be  written  as 


v,  =  X  4,V,  =  ft,  X  «U.M  =  ft,  >),;] 

a-,/=i  *,;=i 

with 

„2  1  _2  (2  +  P)  1 

to  = - ,  to  = - . 

*  Pc  Pls  ‘  P  PcPu, 

The  total  Hamiltonian,  including  the  forcing  current  potential,  is  given  by 
H  =  T  +  VJJ  +  Vn-EJ(2  +  p)ibZ, 

which  yields  (using  Vu  given  in  Eq.  10)  the  equations  of  motion 


(59) 


(60) 


(61) 


C.rr  +  (2  +  p)  C.r  +  2  sin[£  +  YTp  (t  +  Pa  )]  COS(05  +  Tjs  ) 

+  psin[C-i^(0fl+Pj]  =  (2  +  p)4, 


(62) 
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Ms  Vs,rr  +2rls,r  +  2cOS[£  +  ^(0fl  +  7]fl)]sin(f  +  T]s)  +  M s(o]  T}s  =  0, 


(63) 


and 


M 


'),,rr  +  TTplla,:  +  ^sinl^  +  2^(0«  +  *7a  )]«*(&  +  7?s) 


2p 

2+ p 


sin[C  -  +  Va )]  +  Ma(0 1  Tla  =  0. 


(64) 


7.  Solution  for  the  £  Equation  of  Motion  fora  Two-Cell  SQUID  with  /3C  =0 

The  equation  of  motion  for  C  may  be  solved  in  essentially  the  same  manner  as  for  the  single¬ 
cell  SQUID  case.  The  two  sin  terms  are  first  expanded  and  regrouped  to  yield 


M,  £rr  +  (2  +  p)C,T  +  2[sin(£)6!?  +  cos ]  =  (2  +  p)ib, 

with 

=  cos[^(0a  +  I7a)]cos(0,  +  t?s)  +  fcos[^(0fl  +  pj] 

^  =  sin[^(0a  +  J?fl)]cos(05  +  ?7,)  -  f  sin[^(0fl  +  qJJ. 

Writing 

sin(4)a,  +  cos  (£)^  =  [sin(£)cos(7)  +  cos(£)sin(y)]^a  ;  +^2  =  sin(£  +  7?)cf, 

where  _  r1 - 7 

7.  =  tan  (b:  /  aQ,  =  Ja f  , 

the  lowest  order  equation  of  motion  for  C  may  then  be  re-expressed  with  7  =(7f  )  as 

Pc  (C  +  7)>w  +  (C  +  y),T  +  (2cf  /  (2  +  p))  sin(C  +  7)  =4- 

Changing  to  a  new  coordinate 

ir  =  4+r 


(65) 


(66) 


(67) 

(68) 


(70) 


then  yields  an  equation  of  motion  equivalent  to  that  for  a  single  Josephson  junction. 
Setting 


4(2  +  P) 


(71) 


allows  the  equations  of  motion  for  Z  =  0  +  n  /  2  for  the  fc  =0  case  to  be  cast  in  exactly  the 
same  form  as  given  above  in  Eq.  18. 
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Shown  below  are  surface  plots  of  the  critical  current  vs  <fis  and  for  various  values 

of  p ,  and  vs  p  for  the-two  cell  SQUID 


Critical  Current,  p=4 


Figure  12.  Critical  Current  vs.  P  =  l}  Figure  13.  Critical  Current  vs.  {p„,ps.  p  =  4} 


Critical  Current,  (>=03 


Critical  Current,  </>a=2. 


Figure  14.  Critical  Current  vs.  {</>„, <$>„_  p  =  0.3}  Figure  15.  Critical  Current  vs.  {<PS[  P-  <!>„  =  2.0} 
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Critical  Current,  (f)s=2. 


Figure  16.  Critical  Current  vs.  -2.0} 


Pc  =0  Solution  for  the  qs  Equation  of  Motion 

To  derive  the  qs  equation  of  motion,  we  first  re-express  the  cos(£0)and  sin(£0)  terms  as 
cos(£0)  =  cos(^0  -  7)  =  cos(if0)cos(7)  +  sin(i^)sin(y) 

-yjl-a^  sin(ft)0  t)cos(7)  +  (cos(m0  r)  -  d0)sin(y) 

1  -  a0  cos  (cbQ  t) 

sin(£0 )  =  sin(^0  -  y)  =  sin(^0 )  cos(/)  -  cos(^0 )  sin(7)  (72) 

(cos(m0  t)  -  d0)cos(y)  +  yjl-al  sin((y0  T)sin(y) 

1  -  aQ  cos  (ty0  t) 

To  obtain  the  equation  of  motion  to  first  order  in  r\s  we  first  approximate  /]„  — >  (r\a )  ~  0 
to  get 

M,  Vs,n  +  2  Tv,r  +  2  cos(C  +  ^pfl)sin(Ps  +  11 )  +  Msco;  ris  =  0.  (73) 

The  fc  =0  equation  of  motion  to  first  order  in  T]s  may  then  be  written  using  the  results  of  Eq.  53 
as 

+  /(T)[sin(0, )  +  if  cos(0v )]  +  if  /pu  =  0,  (74) 

where 
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and 


y(T)  =  sin(w0T)/sm  +  (cos(w0T)  -  dQ  )/co 


1  -  a0  cos(w0t)  1  -  aa  cos(w0T) 

/sin  =-7^  cos  {Y-k<t>a\ 

/cos=  sin  (y~k<pa\ 


2  +  p 

The  solution  of  this  first-order  ODE  can  be  written  as 

-f  (l/j8Zj+cos(0J)/(T,))dT’  rr  . 

—  — 0  I  Cl 


rj;(T)  =  -e 


fT  ■  IN  /  C/fe+COS^)/^"))^' 

sin(0s)/(T  )e Jo 

•J  () 


r/T'. 


(75) 


(76) 


(77) 


The  integral  in  the  two  exponential  factors  may  be  carried  out  to  yield 


\jvyir 


- sin(y - k(f)a )(ft>0 t) - sjl-aQ  cos (y - k(f)a ) { ln[l - a0 cos(ft)0T)] -in  12} 

sin(y  -  k(j)a  )(1  -  )2  tan-1  [tan((y0T  1 2)  lb] 

a0c°o 


(78) 


The  tan  1  term  above  may  be  re-expressed  as  2  tan  '[tan(ft)0T  12)  l  b]  =00(r)  and  approximated 
by  00(t)  ~  cot)T  .  With  this  approximation  the  solution  can  be  written  as 


rtw  A'. 

[l-a0cos((y(lT')]  [l-a0cos(®0T')]^ 

where 

Q,  =  1  / Pu  +  cos(0, )sin (7  -  k(f)a - 1) /a0, 

and 

^  =  1  ;  cos(0s )  cos(y  -  k<j>a )  =  (2  +  p)  cos(0, )  cos(y  -  k$a ) 

aQih  2  cos  (0J  +  p 


The  final  integration  required  for  the  solution  can  be  carried  out  by  expanding  the 
denominator  in  the  integrand  as  a  series  in  d0  and  can  be  written  as 


(79) 


(80) 


(81) 


(82) 


i!(T)= 


-sin 

d)QdQ(l  —  As) 


f 

,  ..  1,Q,) 

(1  -  a02  )F(t,  X,  ,  Qs )  -  F{ r,  A(  - 1,  Q.s )  1 

J sin 

1  -  aQ  cos(m0r)]  ' 

J  cos 

a0  1  -  aQ  cos((u0t)  1 

where  the  evaluation  of 
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F(t,a,q)=£ 


(83) 


en(T'-T)dz' 

[l  -a0  cos(«0t')]A 

in  terms  of  hypergeometric  functions,  is  presented  in  Appendix  B. 

The  integration  required  for  the  /)"  solution  can  be  carried  out  in  closed  form  if  As  =0  .  This 

will  be  true  if  and  0S  are  near  0.  With  this  assumption  ?7s  can  be  written  as 


n"(r)  =  e0  + 


ec  cos(£90t)  +  es  sin((h0T) 
1  -  a0  cos(g>0t) 


(84) 


where 

e0  =  sin (0, )  sin(y  -  k$a  )a0  /  Q 

Cc  =  -  sin(0s  )[^/l  —  «o  cos(y  -  Wa )  +  sin(7  -  k<t>a  )&]/{(bl+Q,2)  (85) 

es  =  -  sin(0s  )[sin(y  -  )  d>()  -  jl-a2  cos (y  -  k$a  )Q]  /  (©„  +  Q2 ). 


In  the  limit  of  7  — >  0  ,  this  expression  for  q0s  has  the  same  form  as  that  for  a  single-cell 
SQUID  if  the  doubled  cell  area  is  taken  into  account,  i.e., 


(86) 


9.  Pc  =0  Solution  for  the  i)a  Equation  of  Motion 

To  obtain  the  equation  of  motion  to  first  order  in  r\a  we  first  approximate 
cos  (0S  +  r\5)  =  cos  (fs)  =  cos(0s)cos(77s)  -  sin(0s)sin(77s) «  cos(0s)(l  -  (fn; ))  -  sin  (</>s)(r?s)  (87) 
to  get 

Ma  %,rr  +  +  ^sin[£  +  +  7]a)]cOS(# ) 

(87) 

-  l^Sin[C  “2 +  9a  )]  +  %  =  0. 

The  Pc  =0  equation  of  motion  to  first  order  in  i]„  may  then  be  written  as 

nL  +  sin(Co ) H  +  a„rj°a }  +  cos(£o ) {bx  +  bj 7°  }  +  ^  if,  =  0  (88) 

where 
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(89) 


ax  =  cos (k(pa )  cos(0* )  -  cos (K(j)a ) 
bx  =  sin  ma )  cos(0; )  +  sin(X0fl ) 
axx  =-k  sin  (k(j)a )  cos(0s+ )  +  K  sin  (A0fl ) 
=  k  cos  (k(j)a )  cos  (0S+ )  +  A  cos  (K<pa ) 

and  where  we  have  defined 

k  =  — and  A'  =  — — — . 

2  +  p  2  +  p 

With  the  ifa  equation  of  motion  rewritten  as 

<r  +  mt  +  g(r)+(Tt)t  =  0 


a  solution  can  then  be  obtained  as 


0  -\\(2+p)ipu,+f(f)]dr  rt  f  [(2+p)/ft,+/(T")]dT" 

ri°a(T)  =  -e  Jo  -Jo 

•»  0 


*(tW, 


where  the  functions  /(r)  and  g(  T )  are  given  by 


y(r)  =  sin(w0T)/sm  +  (COS(w0T)-fi0)/cos 
1  -  cos(w0t)  1  -  o0  cos(w0t) 

/sm  =  -V1_^o  (fccos(y-fc0Jcos(0j)  + Acos(y  + A0fl)) 

/cos  =  k  sin(y  -  kt )  cos(0; )  +  A-  sin(y  +  A0fl ) 


(90) 


(91) 


(92) 


(93) 


(r)=  sin(w0T)gsm  |  (cos(w0T)-a0)gcos 
1  -  aQ  cos(w0t)  1  —  a0  cos(w0t) 

^sin  =  V1-5 0  (sin(7  -  ^0Jcos(0; )  -  sin(y  +  A0j) 

^cos  =  (cos(r  -  k<Pa )  cos(C )  -  cos(r + Kt )) 


(94) 


The  integral  in  the  two  exponential  factors  may  be  carried  out  as  in  the  fis  case  above  to  yield 


/]»  =-e-n“T 


-L _ f 

c(  rii  t'111  J( 


[l-a0COS(ft)0T')]  “  Jo  [1  -  <2()  cos(cu0T ')] 


—dx\ 


where 


Qfl  =  (2  +  p)  /  fiia  +  {  k  sin(y  -  k(j)a )  cos(0; )  +  A  sin(y  +  A0a )} 


VT^-i) 


(95) 


(96) 
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and 


A  =1- 


k  cos(y  -  k(f)a )  cos(0+ )  +  K  cos(y  +  K(f)a ) 
cos  (0;) 


(97) 


The  final  integration  required  for  the  fi',  solution  can  be  carried  out  by  expanding  the 
denominator  in  the  integrand  as  a  series  in  a0  and  can  be  written  as 


(98) 


<^>sin 

1  „ 

_J_  <^COS 

(l-502)F(T,Afl,Qfl)-F(r,Afl-l,Qfl)l 

®oflo(l _  ^fl) 

l-fioCostthoT)]'  A" 

50 

1  -  d0  cos(th0T) 

W  (t'—t) 


sin(ft)T  ')dr' 


where  the  integrals  J-~c  [1- ficos(®r')]A  and  ■  * 
evaluated  in  Appendix  B. 


wjT'-T)  sin(coT')dz' 

[l-dcos((ur')]A  have  been 
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11.  Appendix  A 

The  integrals  /,  through  7g  are  given  by 

j  =  J_  r2”(cos(x)  — qQ)^  =  Vl-ap  -1 
1  2n^°  l-a0cos(x)  a0 

j  _  1  f2*cos(x)[cos(x)-csin(x)]  ^  _  {2-frg1 +2a0[l-a0(l-a0 +a0b0‘)]}  +  (-1  +  2 a2) 

2k  ^  [1  -  a0  cos(x)]2  2(l-a0)fl02(l  +  fl0)2  2(l-a0)a2(l  +  a0)2.s/l-a02 

_  1  f2;r[cos(x)-csin(x)]  _  a0 
3~2k^  [l-a0 cos(x)]2  X~  (1  -a2)3/2 

j  _  1  r2®  cos(x)3  ^  _  n  [-2o0  +  -  2«q +«q(4  -  5t?0' +2(-l  +  i>01)  + a„(-2+ 6t201)] 

4  _  2tt  Jo  [1-a0cos(jc)]3  ~  +  (l-a0)2(l  +  a0)3 

(—2  +  5^q  +  6aAQ ) 
n  2  Z5J2  J 

(l  ) 
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,  _  1  [cos(x)-csin(x)]2  _  l  +  2a2  +  c2(\-al) 
5_2ttJo  [l-a0  cos(a)]3  A“  2(1-  «q)5/2 
r  =J_f2«  cos(x)  3fl0V 

6  2a-  [1  -  aQ cos(x)]3  '  2(l-fl0)2(l  +  o0)3 


12.  Appendix  B 

The  integral  appearing  in  Eq.  99  may  be  evaluated  by  expanding  the  integrand  denominator  in 
a  series  in  powers  of  a ,  which  yields 

Ax' 

W(r  -r) _ ^ _ 

[l-acos(<yr')]A 

ri>'(,'"{1+^cos«+AW+1|!££ir£a^ ’+w+i)a+2)!^^+...j. 

Each  term  may  now  be  integrated  exactly,  and  the  coefficients  of  each  cos(mwr)  and  sin(mwT)  may 
be  summed  to  yield  a  Fourier  series  with  Fourier  coefficients  proportional  to  a  hypergeometric 
function.  The  final  result  is  given  by 


F(x,X,W)  =  Jr 


[I-Scos(cot)] 


,A  2MV2»  2 


-2T- 


r (n  +  X)f  a  Y  W  cos(n<yr)  +  nwsin(narc) 


Fl(B±T±,!T,n  +  \\a). 


Sr(A)n!UJ  (ncof  +  W1 

Each  term  in  the  summation  increases  the  accuracy  of  the  result  by  roughly  one  decimal  point. 
Integration  by  parts  can  be  employed  to  determine  the  integral 
sm((ox')dx'  _ _ 1 


J1 


[1  -  a  cos(g>t  ')]  cbd(l-X) 


{[1-5 cosCftrr )]*  A - WF( x, X - 1, W )} . 


The  cosine  counterpart  may  be  integrated  directly  to  give 


[1  —  acos((OT  )]  a 
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